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Abstract — The loading of spacecraft propellants is a complex, 
risky operation. Therefore, diagnostic solutions are neces- 
sary to quickly identify when a fault occurs, so that recov- 
ery actions can be taken or an abort procedure can be initi- 
ated. Model-based diagnosis solutions, established using an 
in-depth analysis and understanding of the underlying physi- 
cal processes, offer the advanced capability to quickly detect 
and isolate faults, identify their severity, and predict their ef- 
fects on system performance. We develop a physics-based 
model of a cryogenic propellant loading system, which de- 
scribes the complex dynamics of liquid hydrogen filling from 
a storage tank to an external vehicle tank, as well as the in- 
fluence of different faults on this process. The model takes 
into account the main physical processes such as highly non- 
equilibrium condensation and evaporation of the hydrogen 
vapor, pressurization, and also the dynamics of liquid hydro- 
gen and vapor flows inside the system in the presence of he- 
lium gas. Since the model incorporates multiple faults in the 
system, it provides a suitable framework for model-based di- 
agnostics and prognostics algorithms. Using this model, we 
analyze the effects of faults on the system, derive symbolic 
fault signatures for the purposes of fault isolation, and per- 
form fault identification using a particle filter approach. We 
demonstrate the detection, isolation, and identification of a 
number of faults using simulation-based experiments. 
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1. Introduction 

The loading of cryogenic spacecraft propellants is an inher- 
ently risky and unsafe operation, especially in the case of hy- 
drogen [1-4]. Therefore, diagnostic solutions are necessary 
to quickly identify when a fault occurs, so that recovery ac- 
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tions can be taken or an abort procedure can be initiated be- 
fore system safety is compromised. Model-based diagnosis 
approaches enable quick and robust detection, isolation, and 
identification of faults, because they rely on a detailed model 
of system behavior under nominal and faulty conditions. 

Applying a model-based approach requires an in-depth anal- 
ysis and understanding of the underlying physical processes 
in order to produce an accurate and reliable system model. 
However, cryogenic propellant loading involves complex 
physical processes that are difficult to capture. In this paper, 
we develop a medium-fidelity, lumped-parameter dynamical 
model of propellant loading that takes into consideration a 
variety of complex multi-phase phenomena that govern the 
storage and transfer of cryogenic propellants, yet is simple 
enough to allow for physics analysis and numerical simula- 
tions of real loading systems [5]. We concentrate on a system 
of liquid hydrogen (LH2) filling that is functionally represen- 
tative of the Space Shuttle refueling system. In this system, 
LH2 is stored on the ground in a spherical, insulated, double- 
walled storage tank 1ST), and is transfered to the external ve- 
hicle tank (ET) through a network of pipes and valves. The 
model takes into account the main physical processes such as 
highly non-equilibrium condensation and evaporation of the 
hydrogen vapor, pressurization, and also the dynamics of liq- 
uid hydrogen and vapor flows inside the system in the pres- 
ence of helium gas. Since the model incorporates faults in 
the system, it provides a suitable framework for model-based 
diagnostics and prognostics algorithms. 

We apply a model-based diagnostic approach to the system 
using a combined qualitative -quantitative diagnosis method- 
ology based on the approach of [6]. Deviations in measured 
values from model-predicted values are compared to qualita- 
tive predictions made using the system model for quick fault 
isolation. Fault identification is performed using particle fil- 
ters for joint state-parameter estimation. Using simulation- 
based experiments, we demonstrate the detection, isolation, 
and identification of a number of faults. 

The paper is organized as follows. Section 2 describes the 
propellant loading system and its physics model. Section 3 
overviews the diagnosis approach. Section 4 describes the 
fault detection methodology. Section 5 discusses fault isola- 
tion, and Section 6 develops the fault identification approach. 
Section 7 verifies the approach with a number of simulation- 
based experiments. Section 8 concludes the paper. 
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Figure 1. LH2 propellant loading schematic. 

2. System Modeling 

Model-based diagnosis algorithms utilize a model of the sys- 
tem for the fault detection, isolation, and identification tasks. 
We advocate a physics-based modeling approach, because an 
understanding of the physical processes in both the nominal 
and faulty cases is necessary for successful fault identifica- 
tion. Information on fault severity is necessary in order to 
make appropriate recovery actions in response to a fault, es- 
pecially in the propellant loading domain. In this section, 
we develop the physics model of the LH2 propellant load- 
ing system. Since the focus of the paper is on diagnostics, 
we review the main features of the model, and refer to [4,5] 
for additional details. We first summarize the filling protocol, 
followed by mathematical descriptions of the tank histories. 

Filling Protocol 

The purpose of the LH2 propellant loading system is to move 
LH2 from the ST to the ET. Fig. 1 shows a simplified, but 
functionally equivalent schematic of the system. Initially, 
the ullages of the tanks are at atmospheric pressure due to 
the presence of gaseous hydrogen (GH2). Before filling, the 
tanks are first pressurized. The ST is pressurized to 54.7 psia, 
then 80.7 psia through the use of the vaporizer, which boils 
off LH2 from the ST and returns the GH2 to the ullage of 
the ST. The ET is filled with gaseous helium (GHe) through 
the prepressurization valve, until it reaches 38.7 psia. The 
purpose of pressurization is two-fold. First, it limits potential 
boiling of the propellant by keeping a high vapor pressure in 
the ullage of the tanks. Second, the pressure difference be- 
tween the ST and ET drives propellant from the ST to the ET 
in the absence of a pump. 

Filling progresses in stages with different filling rates, con- 
trolled by the position of the transfer line valve (in reality 
there are a number of valves between the tanks, but in this 
paper we consider a simplified representation consisting of 
a single valve). Slow fill begins first with a low flow rate 
and chilling of the ET. As the liquid drains out of the ST, its 
ullage pressure drops, so the vaporizer constantly maintains 
the ullage pressure to keep LH2 flowing to the ET. The flow 
through the vaporizer valve is modulated based on the error 
between the measured ST ullage pressure and the ST pressure 
set point. As the ET is filled, its ullage volume decreases, and, 
therefore, its ullage pressure increases. The ullage pressure 
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Figure 2. Control Volumes (CV), mass and energy flows in 
an LH2 tank. 


in the ET is maintained using its vent valve, which opens and 
closes to maintain the pressure between 38.7 and 41.7 psia. 

When the ET is 5% full, fast fill begins at a flow rate around 
30 kg/s. When the ET is 72% full, the ullage pressure of the 
ST is reduced to 64.7 psia, reducing the flow rate. When the 
ET is 85% full, the fill rate is reduced further by partially 
closing the transfer line valve. When the ET is 98% full, top- 
ping begins at a lower flow rate. The ET vent valve is also 
opened, reducing the ET ullage pressure to 14.7 psia. Finally, 
at 100% full, topping ends and the tank is then continuously 
replenished until launch to replace the boil off. During re- 
plenish, the fill valve position is modulated to maintain the 
ET level at 100%. 

Tank Modeling 

For each tank, we consider three control volumes: the vapor, 
the liquid, and the vapor film, as shown in Fig. 2. By conven- 
tion, positive mass/energy flows enter the CV, and negative 
flows exit the CV. 

The vapor CV (v) consists of GH2 (subscript v) and GHe 
(subscript g), and is treated as a mixture of ideal gases with 
partial densities p v ( g ) and pressures p v ( g ), as well as a com- 
mon temperature T v , all related to each other by the following 
equations of state: 

Pv — Pv Fly Ty ( 1 ) 

Pg = PgRgTy. (2) 


The liquid CV (l) is where, far from the surface, the tempera- 
ture is equal to 7) and the liquid is treated as incompressible. 
Ti is treated as a constant at 20 K. 

The vapor film CV (/) separates the liquid and gas phases. It 
is treated as saturated hydrogen vapor, whose temperature is 
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equal to that of the liquid/vapor interface, where [1,4,7] 

r ' =Tc te)"- (3) 

with pc = 1.315 MPa, T c = 33.2 K, and n = 5 [7], 

For the liquid CV, the net mass flow is defined by the exter- 
nal mass flow Ji e (i.e., the flow across its boundaries other 
than through the interface), and the interphase (condensation- 
evaporation) flow (defined later): 


where Q v f is the heat flow from the vapor CV, Q fi is the heat 
flow to the liquid CV, and hp, is the enthalpy (heat) of vapor- 
ization. Strictly speaking, hi v depends on the saturated vapor 
temperature, such that it goes to zero when the surface tem- 
perature approaches the critical temperature Tc [9], To take 
this effect into consideration, we use the following simple in- 
terpolation formula for Tf < T c : 

hi v (T f ) = h° lv (|^r) 7 , (ID 


m = J le + Jly (4) 

Similarly, the GH2 and GHe mass conservation for CV (v) is 

TTl-y — Jye. Jlv ( 5 ) 

nig = Jge. (6) 

where and J ge are the external GH2 and GHe flows. For 
a given tank volume V , the vapor volume is fully defined by 
the mass of liquid, since it has a constant density: 

V = V v + V l = V v + m l /p l . (7) 


The energy conservation law for the vapor CV yields: 


d(m v u v + rrigiig) 

dt 


U V Qve Qvf li 


Jivhyf + J ve (h ve + Vy e /2) + J ge (h. ge + v 1 j 2), (8) 


where u v ( g ) = Cv, v ( g )T v are specific internal energies of the 
gases, U v is the net rate of change in internal energy of the 
mixed gas, Q ve is the net external heat flow into the CV (v) 
through the tank walls, Q v f is the heat flow lost through the 
interface, W = — ptdVi/dt is related to the quasi-static power 
due to compression (expansion) of the CV ( v ), h v fJi v is the 
interphase enthalpy flow (the specific enthalpy h v f is equal 
to that of the saturated vapor), h ve and h ge are the specific 
enthalpies of the hydrogen and helium gas entering the CV, 
and v ve and v ge are the velocities of the incoming gases [8], 
Here, the kinetic energies associated with both the GH2 and 
GHe mass flows entering the CV (v) are taken into consider- 
ation, because the corresponding velocities v v ( g ) e are much 
greater than the one related to interphase flow. The specific 
enthalpies h v ( g ) ^vig) T Pv( g )/ Pv(g ) ^v(g) T Ry( g )Ty 

Cp, v (g)T v . Here, c PtV ( g ) = c v , v ( g ) + R v ( g ) and is the specific 
heat at constant pressure. The temperature of the mixed gas 
is then described by 

Ty ( U ‘t! ihyCVyTy TY Ig C \\fj 7); ) , (9) 

m.yCv 

where cy v and Cy, g are the specific heats at constant volume 
for the vapor and gas, and cy is the specific heat at constant 
volume for the mixed gas. 


If the film layer is considered negligibly thin so that one can 
ignore its mass [1], then the energy balance equation for the 
CV (/) can be written as 


Qvf Qfi T Jivhiv — 0 , 


where for liquid hydrogen, Tc = 33 K and h® = u® v = 
4.5 x 10 5 J/kg at p = 1 atm and TJ = 20 K [9]. The heat flow 
terms may be computed based on the liquid, vapor, and film 
temperatures, allowing for to be computed with (10). 

The heat flows Q fi and Q„f are dominated either by conduc- 
tion or convection, depending on the relative temperatures of 
the liquid, vapor, and film [5], If 7/ > T then conduction 
heat transfer defines Q/z, else convection does. If T v >Tf, 
then conduction heat transfer defines Q v f, else convection 
does. In the case of conduction, we use the following approx- 
imation: 

Qf° nd = A f af° nd (T f — T;) (12) 

Q™ f nd = A f a% nd (T v -T f ), (13) 

where Af is the surface area of the interface, and the a terms 
are heat transfer coefficients (see [5]). Typically, conduction 
expressions contain complex integral relations, but here, we 
use this algebraic approximation that has proven adequate for 
our system [5], For convection, we use 

Qfi™ = A f afi nv {T f - T t ) (14) 

Qfi f nv = A f a c ° f nv (T v - T f ). (15) 


Both the liquid and vapor CVs absorb external (e) heat from 
the tank walls. This heat is transferred by means of convec- 
tion so that 

Qv{l)e &v(l)e {Tw T v ppp , (16) 

where A v ^ are the internal tank surfaces in contact with va- 
por (liquid), and a v ^ e are the convection heat transfer coef- 
ficients [5], The wall temperature T w is governed by the heat 
flow passing through the walls from the environment [1-3]. 

The temperature T w of the tank wall, considered uniform, is 
defined by the heat exchange rate with the tank surroundings 
with the effective ambient temperature T a : 

Q w = A w a w (T a - T w ) , (17) 

where A w is the external tank surface area. The wall temper- 
ature is governed by the tank energy conservation: 

TTl W Gyj Tyj : Q u ! Qle Q ve- (18) 


( 10 ) 
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Here, the heat transfer coefficients describe natural convec- 
tion inside and outside the tank walls [1,10]. 



The Storage Tank 

The above equations apply equally well to both the ST and 
the ET. We denote variables of the ST with a 1 subscript, and 
variables of the ET with a 2 subscript. For the ST, there is 
no GHe in the ullage volume, so p gl = (> g \ = 0. The exter- 
nal mass flow for the liquid consists of vaporizer flow J vap , 
transfer line flow J tr , and leak flow J j ea k:\ ■ 

Jlel — Jvap Jtr Jl.leakl * (19) 

J tr is described by 

Jtr — Xt r G t.rtXtrs/ P\ P 2 t ( 20 ) 


evaporation as the ET walls are being initially chilled down 
during the beginning of the slow fill stage. 


In the vapor CV, the external flows are given by 


Jve 2 Jv,vent,2 Jv,leak,2 

Jge 2 Jg.pp Jv, vent, 2 Jg.leak. 2i 

where J g . vp is the prepressurization flow of GHe, and 


(26) 

(27) 


Jv(g) ,vent,2 ent ,2 rr V ent .2 


Avent ,2pv(g) \J 't(.Pv2~\~Pg2 Patm ) 


r \J Krent,2(Pv2+Pg2) 


(28) 


where X tr is the valve position fan input), at r is the flow 
coefficient, pi = p v \ + pighn is the total ST pressure where 
hn is the height of the liquid, p2 = p v 2 +Pg2 + Pighi 2 is the 
total ET pressure, g is the acceleration due to gravity and a tr 
is a (dimensionless) multiplicative factor describing a transfer 
line blockage fault. 

In the vaporizer, a certain amount of LH2 is evaporated and 
returned to the ST, thus controlling its ullage pressure. We 
assume that all liquid flow through the vaporizer is converted 
instantaneously to GH2. The flow is given by 


Jvap — ''•vap® vap^vapyj Pv 1 Pvap 

where X vap is the vaporizer valve position, a vap is the flow 
coefficient, p vap is the pressure in the vaporizer (close to at- 
mospheric pressure), and cr tr is a multiplicative factor de- 
scribing the vaporizer valve blockage fault. The vaporizer 
valve position is controlled by 

X vap = min ^1, max ^0, W Pv1 ^ Pvl ^ ^ , (22) 

where p* vl is the desired ST ullage pressure. 

The liquid leak flow is given by 

Jl,leak,l — Aii ea k,lCXpi ea kpy/pi Patm? (23) 

where Ai t i ea kp is the the area of the leak hole, api ea k,i is the 
leak hole flow coefficient, and p a tm is atmospheric pressure. 

In the vapor CV, J„ e i = —J v ,ieak, 1 , where 

Jy,leak,l = A v ,i ea k pCH-v .leak V Pvl Patm ■ (24) 

Heat leaks constitute the external heat flow term, i.e., Q ve i = 

Qleak.l • 

The External Tank 

For the ET, the external liquid flow is given by 

Jle 2 ~ Jtr Jl,leak.2 Jboili (25) 

where Jpi ea k 2 is described similar to that for the ST, and 
Jboii = Qie 2 /hi v (T W 2 — T 12 ) is responsible for intense LH2 


Here, the dimensionless flow coefficient K (the loss factor) 
can be found in Schmidt et. al. [10] (see Tables 7-2 and 7-3 
therein); a dimensionless relative valve position assumes val- 
ues between Afc = 1 (fully open) and A/ c = 0 (fully closed) 
[10], During filling, the valve opens when the pressure ex- 
ceeds 41.7 psia and closes when it falls below 38.7 psia. The 
vapor/gas and heat leaks may be defined as with the ST. 

Nominal Dynamics 

Fig. 3 summarizes the major results of the simulation of a 
nominal loading regime, which are based on the parameters, 
initial conditions, and filling protocol that are typical for LH2 
loading systems. These values are provided in [5]. 

It can be seen that the LH2 level in the ST drops monoton- 
ically (Fig. 3a) as the level in the ET rises (Fig. 3d). The 
pressure p v i in the ST (Fig. 3b) is determined by the load- 
ing dynamics and controlled by the vaporizer. Once achieved 
during slow fill, the pressure in the ST is maintained at ap- 
proximately 80.7 psia up to the start of the reduced-pressure 
fast fill (see Figs. 3b and 3j), at which point it is maintained 
at 64.7 psia. Meanwhile, the ET ullage pressure (Fig. 3e) is 
oscillating due to the cycling of the vent valve that maintains 
the pressure between the lower and upper thresholds of 38.7 
and 41.7 psia. The fluctuations in the ET ullage temperature 
(Fig. 3f) as well as in the mass flow rates (Figs. 3j and 3h) are 
driven by the ET pressure oscillations. 

The LH2 partial pressure in the ET rises due to the contin- 
uing hydrogen supply, while the GHe partial pressure drops 
because the helium is being permanently removed through 
the vent valve (Fig. 3i). In this case, due to the condensation 
blocking effect [4], the flow of the condensed vapor in the ET 
(Fig. 3h) is much smaller than that in the ST (Fig. 3g), be- 
cause the vapor pressure is being maintained approximately 
equal to the equilibrium pressure of the condensed vapor at 
the temperature of LH2. The ullage temperature T v 2 in the 
ET (Fig. 3f) increases initially due to the introduction of the 
GHe during the pressurization stage, then drops from the ini- 
tial high value due to venting and near-wall boiling that gen- 
erates relatively cold GH2 during filling. The liquid surface 
temperature Tf 1 in the ST increases (Fig. 3c) due to the vapor 
condensation at the interface. Simultaneously, the ST ullage 
temperature T v 1 increases, mainly because the relatively hot 
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Figure 3. Nominal regime of LH2 propellant loading. 


GH2 is supplied by the vaporizer as loading occurs. As a 
result, the ullage temperature approaches the temperature of 
LH2 saturated vapor at a pressure close to the final ST ullage 
pressure of approximately 5 atm (Fig. 3b). 

3. Diagnosis Approach 

We apply a model-based diagnosis approach using the 
physics model of the LH2 system. In this paper, we con- 
sider the problem of single fault diagnosis. The system may 
be described in the following general form: 

*(t) = f(t,x(t),0(t),u(t),v(t)) 
y(t) = h(t,x(t),0(t),u(t),n(t)), 


where x(t) G R” x is the state vector, Q(t) € R " 8 is the 
parameter vector, u(i) € R"“ is the input vector, v(<) € R" 1 ’ 
is the process noise vector, f is the state equation, y(i) € R" y 
is the output vector, n(f) € R n " is the measurement noise 
vector, and h is the output equation. 

Measurements are time-varying signals of y(t) obtained from 
the system sensors. In the LH2 system, we consider the fol- 
lowing measurements for diagnosis: Vn, pi, T„i, V 12 , P 2 , 
Jtn ^vap-> an d A ven t, 2 ■ 

We consider single, abrupt faults, modeled as unexpected step 
changes in system parameter values. We name faults by the 
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Figure 4. Diagnosis architecture. 


associated parameter and the direction of change, i.e., 0 de- 
notes a fault defined as an increase in the value of parameter 9 , 
and 9~ denotes a fault defined as a decrease in the parameter 
value. For the LH2 system, we consider liquid and gas leaks, 
heat leaks, and valve clogging. The liquid and vapor leaks are 
defined by the equivalent leak areas, nominally zero, so faults 
are defined as increases in these areas, denoted by Af, , 1 , 
A + , , , , Af, , „, and A + , , „. The heat leaks are defined 
by the heat leak rate, nominally zero, so faults are increases 
in these values, denoted by Qf eak x and Qf eak 2 - The valve 
clogging faults are described by the a parameters, which are 
nominally 1, so faults are decreases in these values (0 at a 
minimum), denoted by a^ r , a~ ap , and cr“ ent 2 - 


residual signal, y, r (t), is estimated over a small window W 2 : 

= w 2 E r ^- 

z i=t-w 2 + 1 

The variance of the nominal residual signal, is com- 

puted using a large window W\ preceding W 2 , by a buffer 
Wdeiay ensuring that W\ does not contain any samples after 
fault occurrence. The variance is computed using 

t-w 2 -w delay 

of (t) = TTT- ( r W “ 

1 i^t-Wi-W^y-Wt+l 

where 


The diagnosis architecture is shown in Fig. 4. The system 
receives inputs u(f) and produces outputs y (t). The physics 
simulation runs simultaneously, producing predicted outputs 
y(t), given the inputs u(f). Using statistical methods, the 
fault detection module decides when a measurement has de- 
viated from its nominal value, triggering fault isolation. Mea- 
surement deviations are then used to quickly isolate faults F. 
Fault identification computes, for each fault / £ F, the value 
of the fault parameter that best fits the outputs of the system, 
and the candidate with the lowest output error is selected as 
the best candidate. 

4. Fault Detection 

In model-based fault detection, a model of the system pro- 
vides reference outputs representing nominal system behav- 
ior. For each sensor output y(t), we define the residual as 
r(f) = y(t) — y(t), where y(t) is the model-predicted out- 
put signal. Statistically significant deviations of the actual 
system outputs from the model-predicted outputs imply the 
presence of a fault. If the model is accurate, then fault detec- 
tion thresholds can be small and faults detected quickly. The 
thresholds are dynamic in that they are defined with respect 
to nominal behaviors as a function of time. This is favorable 
to the current practice with such systems, where thresholds 
are static, preventing the detection of subtle deviations from 
nominal behavior that indicate faults. 

We use the Z-test for robust fault detection using a set of slid- 
ing windows, as described in [ 1 1, 12], The current mean of a 


/4W 


1 

IUi 


t—W2~W delay 

r (*)- 

i=t—W 2 — Wdeiay ~Wl+l 


A given confidence level determines the bounds z~ < 0 and 
z + > 0 for a two-sided Z-test. The fault detection thresholds, 
E~ (f ) and sf(t), are dynamically computed using 


£ r ( t ) 
£f(t) 


' - CTrjt) 
Z VW 2 

CT r (t ) 


-E 
+ E, 


where E is a modeling error term. A fault is detected if /x r (<) 
lies outside of the thresholds at time t. If (£)> a 

- symbol for the measurement is used by the fault isolation 
module, and if y, r (t) > sf{t), a + symbol is used. Generally, 
the parameters W\, W 2 , Wdeiay , the z bounds, and E must be 
tuned to optimize performance to minimize both false alarms 
and missed detections. 


5. Fault Isolation 

We utilize a qualitative diagnosis methodology that isolates 
faults based on the transients they cause in system behav- 
ior, manifesting as deviations in observed measurement val- 
ues from nominal measurement values [6], The transients are 
abstracted using qualitative + (increase), - (decrease), and 0 
(no change) values to form fault signatures. Fault signatures 
represent these measurement deviations from nominal behav- 
ior as the immediate (discontinuous) change in magnitude, 
and the first nonzero derivative change. 


6 








The fault signatures can be derived automatically from a 
graph-based representation of the system model, known as a 
temporal causal graph (TCG) [6]. A TCG captures dynamic 
system variables as nodes in a graph, and the qualitative re- 
lations between them as edges. These edges are labeled with 
dt, representing integration, +1, denoting a proportionality, 
—1, denoting an inverse proportionality, and system parame- 
ters. Faults that are modeled as parameter changes appear on 
edges, allowing the qualitative effects of parameter changes 
to be propagated over the system variables. Propagation to 
measured variables reveals the qualitative effects of faults on 
measurements. 

A partial TCG for the LH2 system, showing some key vari- 
ables of the ST, is shown in Fig. 5. Variables within dashed 
boxes are measured. Here, we show only the cr vap fault. A 
forward propagation algorithm may be used to derive the fault 
signatures. Details may be found in [6], so here, we use the 
vaporizer valve clogging fault as an example to illustrate the 
general procedure. According to the TCG, a decrease in a vap 
will lead to a decrease in J vap , and, subsequently, a decrease 
in and, due to the —1 label, an increase in Jn . The de- 
crease in J v i will lead to a first-order decrease (due to the 
dt label) in m v i and, subsequently, p v i, which is a measured 
variable. This change then propagates as a decrease in J tr , 
and then towards ET variables. The increase in Jn will lead 
to a first-order increase in mu and Vn, which is measured, 
leading to an increase in J tr . This conflicts with the previ- 
ously predicted decrease from the other path with the same 
derivative order, resulting in an ambiguity denoted with a * . 
We can use the simulation and knowledge of the system dy- 
namics to resolve these ambiguities. In this case, we know 
that the hydrostatic pressure due to the height of the liquid 
is negligible compared to p v i due to the very low density of 
LH2. Therefore, we know that the decrease effect will domi- 
nate. The propagation continues to all system variables. 

Fault signatures for the LH2 system are shown in Table 1. 
Due to the dynamics of the system, all faults appear as smooth 
changes in the measured values, so only the first change 
presents useful diagnostic information, hence, we show only 
the first symbol derived from forward propagation. Recall 
that the derived symbols represent deviations from nominal 
behavior, so, for example, the + symbol for Vn caused by 
leak i does not necessarily mean that Vn increases, rather, 
it increases with respect to its nominal value, i.e., it decreases 
at a slower rate. The remaining * symbols are for those in 
which the qualitative effect may be different depending on the 
state of the system. Note also that this is a hybrid system, i.e., 
it has continuous dynamics mixed with discrete dynamics due 
to the valves, and we show only the signatures for the mode 
where the vent valve is closed. When the vent valve opens, 
some of the signatures will flip. A systematic framework for 
dealing with this issue is discussed in detail in [13]. 

Using the signatures, we can perform diagnosability analy- 
sis to determine the effectiveness of the isolation step. We 
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Figure 5. Partial TCG of the LH2 system. 


Table 1. Fault Signatures for the LH2 System 


Fault 

Vn 

Pi 

T v i 

Vn 

P2 

T v 2 

Jtr 

A vap 

^vent, 2 

4 + 

^l,leak, 1 

- 

- 

* 

- 

- 

* 

- 

+ 

- 

4 + 

^v,leak ,1 

+ 

- 

- 

~ 

~ 

* 

- 

+ 

- 

Q leak ,1 

- 

+ 

+ 

+ 

+ 

■ k 

+ 

- 

+ 

4 + 

leak, 2 

- 

- 

•k 

~ 

~ 

■ k 

+ 

+ 

- 

A + 

^v,leak, 2 

- 

- 

•k 

+ 

- 

- 

+ 

+ 

- 

Q leak ,2 

+ 

+ 

■ k 

~ 

+ 

+ 

- 

- 

+ 

a tr 

+ 

+ 

■ k 

~ 

~ 

~ 

~ 

- 

- 

®vap 

+ 

~ 

~ 

~ 

~ 

■ k 

~ 

+ 

- 

G vent ,2 

+ 

+ 

* 

~ 

+ 

+ 

~ 

- 

+ 


can see that in most cases, the set of signatures produced dis- 
tinguishes most faults. One exception is the pair A' v leak 1 
and crj ap , which have all the same signatures, therefore, fault 
identification will have to distinguish between them, since 
quantitatively, they should have different effects on the mea- 
sured variables. If we measure J vap , then the faults can 
be distinguished, because aj ap would produce a decrease 
in Jvap , whereas A+ leak 1 would produce an increase. The 
faults Qn ak o and a~ ent 2 are similarly indistinguishable. 

We can also use the signatures to perform measurement se- 
lection. The measurement of \ vap always has the same sig- 
nature as pi, therefore one of these can be eliminated with- 
out loss of diagnosability. However, both quantities are rela- 
tively easy to measure, so there is value in keeping them both. 
Xvap is typically more sensitive to faults than p \ , since p\ is 
controlled by the vaporizer, which may mask faults. Simi- 
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larly, both p 2 and A ven tp are not needed. Changes in P 2 are 
reflected in changes in the switching frequency of the vent 
valve, and therefore provide the same information, although 
changes in the vent valve position can be detected much more 
easily, again because the control of p 2 leads to fault masking 
in that measurement. It is also beneficial to keep as many 
measurements as possible, because this gives more informa- 
tion for fault identification. Further, the qualitative fault sig- 
natures say nothing about fault magnitude. Some faults may 
not produce large enough changes in some variables for the 
changes to be detected, therefore, having more measurements 
helps to alleviate this problem. 

6. Fault Identification 

Fault identification is initiated immediately after the initial set 
of fault candidates is produced after fault detection. Each can- 
didate has its own identification module that updates its esti- 
mate at every time step. Identification is performed using the 
particle filtering algorithm for joint state-parameter estima- 
tion [14, 15], which has seen previous application in model- 
based diagnosis and prognosis algorithms [16-19]. The mag- 
nitude of the fault parameter 0 is estimated along with the 
system state. 

The identification module must compute p(x.k, 0fc|yo:fc)- A 
general solution to this problem is the particle filter, which 
may be directly applied to nonlinear systems with non- 
Gaussian noise terms. Particle filters offer approximate (sub- 
optimal) solutions to the state estimation problem for systems 
where optimal solutions are unavailable or intractable [14, 
15]. In particle filters, the state distribution is approximated 
by a set of discrete weighted samples, called particles. As 
the number of particles is increased, performance increases 
and the optimal solution is approached. Due to the highly 
nonlinear dynamics of the LH2 loading system, particle fil- 
ters are favored over other estimation approaches such as the 
extended Kalman filter. 

The particle approximation to the state distribution is given 
by 

{( x fc>0fc)>Wfe}£i> 

where N denotes the number of particles, and for particle i, 
x], denotes the state vector estimate, 0 ). denotes the parameter 
vector estimate, and ufi denotes the weight. The posterior 
density is approximated by 

N 

p(x fe ,0 fe |y 0: fc) « Wk&(xifii)(dx-kd9k ), 

i = 1 

where <5 (x i ni\(dxf c dOf c ) denotes the Dirac delta function lo- 
cated at (x|,0* fe ). 

We use the sampling importance resampling (SIR) particle 
filter with systematic resampling [14,20], The pseudocode 
for a step of the filter is shown as Algorithm 1 . Each particle 


Algorithm 1 SIR Filter 

Inputs: {(x^_ 1 , 0j._ 1 ), "4-tJiLi, u k~i:k, y/c 
Outputs: «,*}£! 

for i = 1 to TV do 

ei~p(efc|e , fc - 1 ) 

x* ~p(x fc |x i fc _ 1 ,e|_ 1 ,u fe _ 1 ) 

w i <- p(yk\x-i,0i,uk) 

end for 

N 

w^Y. w i 

i= 1 

for i = 1 to N do 

w k W H W 

end for 

{( x fc> e fc)> w fc}£i Resample({(x*,0j.),w*}y x ) 


is propagated forward to time k by sampling new parameter 
and state values. The particle weight is assigned using y^. 
The weights are then normalized, and then the particles are 
resampled (see [14]). 

Note that the parameters 6 *. evolve by some unknown process 
that is independent of the state x^. The particle filter algo- 
rithm requires some type of evolution to the parameters. We 
use a random walk, i.e., for parameter 0, Ok = Ok - i + £fc-i, 
where Gc-i is Gaussian noise. The particles generated with 
parameter values closest to the true values should match the 
outputs better, and, therefore, be assigned higher weight, thus 
allowing the particle filter to converge to the true values. The 
selected variance of the random walk noise affects both the 
rate of this convergence and the estimation performance after 
convergence. Using the simulation, we can determine appro- 
priate values for the random walk variances. 

Under the single fault assumption, we can run a set of par- 
allel particle filters, one for each consistent fault candidate. 
This reduces the dimensionality of the estimation task over 
combined estimation of all fault parameters, and allows the 
particle filters to be much more efficient. The particle filters 
are initiated at the point of fault detection, using the model- 
estimated state at that time point for initialization of the states. 
The fault parameter estimate starts at its value during nomi- 
nal operation (e.g., 0 for leak areas). When the fault isolation 
module reduces the set of candidates F, the identification task 
continues only for those faults remaining in F. 

Because diagnosability may be limited, fault identification 
must also be used to help refine fault candidates. The par- 
ticle filter for the true fault candidate should estimate the cor- 
rect fault parameter and track the faulty outputs with low er- 
ror, whereas the particle filters for the incorrect faults will not 
track and result in large error. We compute the mean squared 
output error for each candidate from the point of fault detec- 
tion to the present time. The candidate with the lowest output 
error e is considered to be the true candidate. 



7. Results 


We illustrate the diagnosis process with the A+, , , fault 
injected at 2600 s with magnitude 3 x 10 -3 m 2 . Note that 
from diagnosability analysis, the observed signatures will be 
consistent with both the ST vapor leak fault A plk 1 and the 
vaporizer valve clogging fault a~ ap , so the fault identification 
stage will have to resolve the ambiguity. Relevant measure- 
ments are shown in Fig. 6. The diagnoser, running at a sample 
time of once per second, detects the fault at 2605 s due to an 
observed decrease in p\. The initial list of consistent faults is 

leak, 1' -^v, leak, 2' A + Zeafc,l > leak, 2’ a vap}' At 2608 S, a 

decrease in J tr is detected, followed by an increase in X vap 
at 2612 s and a decrease in T v i at 2690 s. The candidate 
set remains the same. At 2703 s, a decrease in V 12 is de- 
tected, which rules out the ET vapor/gas leak fault A+ leak 2 . 
At this point, the output error of the At, , 1 fault is an or- 
der of magnitude less than for all other remaining candidates. 
At 2806 s, a detected increase in Vn rules out Af leak 1 and 
Kieak^ resulting in the candidate list {A+ leak l ,a~ ap }. By 
this time, the particle filter for At, k 1 has also converged on 
the correct value of the fault magnitude, as shown in Fig. 7. 
The fault was detected within 5 seconds of its occurrence, and 
identification confirmed the true fault with correct magnitude 
within 200 s. It is important to quickly discriminate between a 
vapor leak in the ST and a clogging of the vaporizer, because 
the former requires an abort, whereas the system can continue 
fueling in a safe manner with the latter. As discussed in Sec- 
tion 5, measuring J vap would allow a faster discrimination of 
these two faults. 

Diagnosis results over the complete set of faults are shown in 
Table 2. Here, At ( i denotes the time to detect the fault. At, 
denotes the time to isolate, taken as the last time at which 
the candidate set is reduced. Fid denotes the output of the 
fault identification module at the end of the scenario, and f* 
denotes the final output of the diagnoser, i.e., the identified 
fault with the lowest error. For the fault identification stage, 
50 particles were used per particle filter, which seemed to 
offer a reasonable trade off between computation time and 
identification accuracy. 

Overall, the results are fairly good. In all cases, the correct 
fault was identified with sufficient accuracy. In some cases, 
the effects of the faults are very subtle at first, and it takes 
some time before the changes they produce can be distin- 
guished from the sensor noise. For the vent valve clogging 
fault, the effects are visible only when the valve is open, so 
if the fault first appears when the valve is closed, it will take 
time for the fault to be detected. The erjj. fault was detected 
immediately because the change in J tr was significant. With- 
out this measurement, it would take some time for the fault to 
be visible in the liquid volumes of the tanks. 

For some faults, the operation of the vaporizer masks the 
change in ST pressure, therefore, it is important to monitor 
the vaporizer valve position in order to obtain a more sensitive 



2000 2500 3000 3500 4000 



2000 2500 3000 3500 4000 

t (s) 

Figure 6. Predicted and observed outputs for the A+ , . 1 
fault injected at 2600 s with magnitude 3 x 10 -3 m 2 . 



2600 2800 3000 3200 3400 3600 



2600 2650 2700 2750 2800 2850 2900 2950 3000 

t (s) 

Figure 7. Estimated fault magnitudes and output error for the 
A+ leak j fault injected at 2600 s with magnitude 3 x 10” 3 m 2 . 
The A pl k | estimator converges to the true value by 2800 s. 
Thecr^op estimator converges by 2800 s also, but it cannot 
track the outputs, as shown by the corresponding error e& vap , 
which is greater than ; i at 2800 s. 


9 



Table 2. Diagnosis Results 


Fault 

Magnitude 

73 

-to 

<1 

At; 

Eid 

f* 

4 + 

^l,leak, 1 

1 X icr 3 m 2 

17 

113 

A Lleak. 1 = 9.98 x 1 0-4,e = 97.1 

A t,eak. 1 = 9 ' 98 >< 10 ~ 4 

4 + 

v , leak, 1 

3 x icr 3 m 2 

5 

206 

<ieafc, i = 3 - 04 x 10— 4 ,e = 27.1 
Gvap = 0.35, e = 7.89 x 10 5 

A t,leak,l = 3 -° 4 >< iO " 4 

Qleak, 1 

3 x 10 4 W 

111 

151 

Qteak. i = 2.95 Xl 0 4 ,e = 9.44 

Qteak. 1 = 2 ' 95 X 10 4 

4 + 

^l,leak, 2 

1 x 10“ 3 m 2 

28 

131 

A t.leak. 2 = E10xl0~ 3 ,e=121 

A t.leak. 2 = 1-10 X 10~ 3 

A + 

,leak,2 

1 x icr 4 m 2 

39 

138 

Kleak.i = 1 - 30 X 10- 4 ,e = 409 

A t,leak, 2 = 1-07 X IQ " 4 


Kleak, 2 - !.°7 x !°-V = 7.50 
Kleak, 1 = 6.78 Xl0- 5 ,e = 463 
^m = U6x 10 - 4 ,e = 2.S6 x 105 
ovT.,,, = 0.95, e = 410 


Qteak, 2 

3 x 10 4 W 

18 

131 

QL M = 4 - 15xlo3 ’ e = 

<?L fe .2 = 3 - 19 xl0 4 ,e = 

= 1930 
= 8.97 

Qteak , 2 = 319 X 10 4 

a fr 

0.5 

0 

13 

o-- = 0.50, e = 56.8 


af r = 0.50 

&vap 

0.5 

7 

71 

A t,leak,l = 2 - 14 X 10- 3 , 
crjap = 0.50, e = 8.06 

e = 3.47 x 10 5 

&vap = 0-50 

®vent,l 

0.5 

141 

152 

QL fe ,2 = 3 - 21xlo3 ’ e = 

^■ e „t.l = °- 50 - e = 9 - 55 

= 105 

i = 0-50 


detection. Changes in the vent valve switching frequency are 
also easy to observe, but one must wait for the upper pressure 
threshold to be reached before it can be determined whether 
it closes early or late. We also observed that faults in the ET 
have very small effects on measured variables in the ST, re- 
sulting in the large size of for A+ leak 2 . In such cases, 
fault identification becomes even more crucial to resolving 
ambiguities. 

8. Conclusions 

In this paper, by applying our previously developed physics- 
based model of a propellant loading system [5], we analyzed 
the effects of faults, and applied a model-based diagnosis ap- 
proach to fault detection, isolation, and identification. The 
detection stage compares observed and model-predicted out- 
puts to detect faults. The isolation stage compares model- 
predicted fault transients to observed measurement deviations 
to quickly isolate faults and reduce the complexity of the 
identification stage. The fault identification stage uses par- 
ticle filters to estimate the values of fault parameters and re- 
solve isolation ambiguities. 

The simulation results showed that with such a model-based 
diagnosis approach, faults can be quickly detected, isolated, 
and identified, and this knowledge can be used to determine 
if the system can continue loading safely or if an abort is re- 
quired. We considered only single faults here, and in the fu- 
ture we will extend the framework to multiple faults. Sensor 
faults were not considered here, but can be easily incorpo- 
rated into a model-based framework by including models of 
the sensors [12], Using the physics model for prognostics and 
loading optimization is also of interest. 

We would like to validate the approach by using available his- 
torical data of faulty situations. Such data can be difficult to 


find, and, further, there is only a limited capability to inject 
faults into the system for the purpose of diagnosis algorithm 
validation, due to the high cost of a fueling operation and the 
explicit danger that fault injection presents in such a system. 
Therefore, an accurate simulation with which to validate di- 
agnosis algorithms, such as that used here, is very valuable. 
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